function D = shockdb(varargin)
% shockdb  Create model-specific database with random shocks.
%
%
% Syntax
% =======
%
% Input arguments marked with a `~` (tilde) sign may be omitted.
%
%     D = shockdb(M,~D,Range,~NDraw,...)
%
%
% Input arguments
% ================
%
% * `M` [ model ] - Model object.
%
% * `D` [ struct | empty ] - Input database to which shock time series will
% be added; if omitted or empty, a new database will be created; if `D`
% already contains shock time series, the data generated by `shockdb` will
% be added up with the existing data.
%
% * `Range` [ numeric ] - Date range on which the shock time series will be
% generated and returned; if `D` already contains shock time series
% going before or after `Range`, these will be clipped down to `Range` in
% the output database.
%
% * `NDraw` [ numeric ] - Number of draws (i.e. columns) generated for each
% shock; if omitted, the number of draws is equal to the number of
% alternative parameterizations in the model `M`, or to the number of
% columns in shock series existing in the input database `D`.
%
%
% Output arguments
% =================
%
% * `D` [ struct ] - Database with shock time series added.
%
%
% Options
% =======
%
% * `'randFunc='` [ `@lhsnorm` | `@randn` | *`@zeros`* ] - Function used to
% generate random draws for new shock time series; if `@zeros`, the new
% shocks will simply be filled with zeros; the random numbers will be
% adjusted by the respective covariance matrix implied by the current model
% parameterization.%
%
% Description
% ============
%
%
% Example
% ========
%

% -IRIS Toolbox.
% -Copyright (c) 2007-2015 IRIS Solutions Team.

[This,D,Range,NDraw,varargin] = ...
    irisinp.parser.parse('model.shockdb',varargin{:});
opt = passvalopt('model.shockdb',varargin{:});

%--------------------------------------------------------------------------

ne = sum(This.nametype==3);
nPer = length(Range);
nAlt = size(This.Assign,3);
eName = This.name(This.nametype==3);
eLabel = mylabelorname(This);
eLabel = eLabel(This.nametype==3);

if isempty(D)
    E = zeros(ne,nPer);
else
    E = datarequest('e',This,D,Range);
end
nShock = size(E,3);

doChkSize();

nLoop = max([nAlt,nShock,NDraw]);
if nShock==1 && nLoop>1
    E = repmat(E,1,1,nLoop);
end

randFunc = func2str(opt.randfunc);
switch randFunc
    case 'lhsnorm'
        S = lhsnorm(sparse(1,ne*nPer),speye(ne*nPer),nLoop);
    otherwise
        S = opt.randfunc(nLoop,ne*nPer);
end

for iLoop = 1 : nLoop
    if iLoop <= nAlt
        Omg = covfun.stdcorr2cov(This.stdcorr(1,:,iLoop),ne);
        F = covfun.factorise(Omg);
    end
    iS = S(iLoop,:);
    iS = reshape(iS,ne,nPer);
    E(:,:,iLoop) = E(:,:,iLoop) + F*iS;
end

% `E` is ne-by-nPer-by-nLoop, permute to nPer-by-nLoop-by-ne.
E = permute(E,[2,3,1]);

temp = tseries();
for i = 1 : ne
    name = eName{i};
    D.(name) = replace(temp,E(:,:,i),Range(1),eLabel{i});
end


% Nested functions...


%--------------------------------------------------------------------------


    function doChkSize()
        if nAlt>1 && NDraw>1 && nAlt~=NDraw
            utils.error('model:shockdb', ...
                ['Input argument NDraw is not compatible with the number ', ...
                'of alternative parameterizations in the model object.']);
        end
        
        if nShock>1 && NDraw>1 && nShock~=NDraw
            utils.error('model:shockdb', ...
                ['Input argument NDraw is not compatible with the number ', ...
                'of alternative data sets in the input database.']);
        end
        
        if nShock>1 && nAlt>1 && nAlt~=nShock
            utils.error('model:shockdb', ...
                ['The number of alternative data sets in the input database ', ...
                'is not compatible with the number ', ...
                'of alternative parameterizations in the model object.']);
        end
    end % doChkSize()
end
